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Abstract 

Spatiotemporal simulation of minimum and maximum temperature is a fundamen- 
tal requirement for climate impact studies and hydrological or agricultural models. 
Particularly over regions with variable orography these simulations are difficult to pro- 
duce due to terrain driven nonstationarity. We develop a bivariate stochastic model for 
the spatiotemporal field of minimum and maximum temperature. The proposed frame- 
work splits the bivariate field into two components of "local climate" and "weather". 
The local climate component is a linear model with spatially varying process coeffi- 
cients capturing the annual cycle and yielding local climate estimates at all locations, 
not only those within the observation network. The weather component spatially cor- 
relates the bivariate simulations, whose matrix valued covariance function we estimate 
using a nonparametric kernel smoother that retains nonnegative definiteness and allows 
for substantial nonstationarity across the simulation domain. The statistical model is 
augmented with a spatially varying nugget effect to allow for locally varying small scale 
variability. Our model is applied to a daily temperature dataset covering the complex 
terrain of Colorado, USA, and successfully accommodates substantial temporally vary- 
ing nonstationarity in both the direct-covariance and cross-covariance functions. 

Keywords: complex terrain, Gaussian process, nonstationary, minimum 
temperature, maximum temperature, multivariate covariance, simulation, 
stochastic weather generator 
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1 Introduction 



Stochastic simulation of physical variables such as minimum or maximum temperature, pre- 
cipitation amount and solar radiation are often required as inputs to physical models over 
varying types of topography. Over plains regions, agricultural and crop models require daily 
minimum and maximum temperature simulations at locations that typically do not have 
direct observations. In mountainous regions, hydrological models require stochastic weather 
realizations for runoff, snowmelt and watershed modeling, as well as water resourc e planning 



and climate impact assessment uKustas et al 



1994 : 



Semenov and Barrow 



1993). 



Stochastic weather generators (SWGs) are one approach to producing simulations of daily 



weather; they 



1996 



are simply prob a 



obse rvations ( 



Wilks and Wilbv 



Racsko et al 



1991 



jility models whose simulations are statistically similar to 



Richardson 



Rajaaopalan and Lall 



1999). SWGs can loosely be categorized into model-based 



198ll ) and empirical approaches (e.g. 



Lall and Sharma 



1999). Often these weather generators produce simulations 
only at locations with observational data, but modern physical models require gridded daily 
weather. Hence, recent research has been directed toward gen erating spatially consistent 



SWGs that are available at and between observation locations ( Wilks 



1999 



Kleiber et al. 



20121 ). Herein we focus on a model-based approach to minimum and maximum temperature 



simulation over a mix of complex terrain and relatively homogeneous terrain simultaneously. 

Spatially consistent simulation over most agricultural regions can be accommodated using 
isotropic or stationary models that are appropriate for regions with relatively constant or 
slowly changing topography. Domains with highly variable terrain, in particular mountain- 
ous domains, are challenging for the majority of univariate spatial models due to substantial 
nonstationarity of physical processes in these areas. Weather over complex terrain is highly 
variable due to topography; for example, at high elevations in the northern hemisphere, north 
facing slopes tend to be cooler than lower elevations and south facing slopes and valleys can 
create their own micro-climate relative to the surrounding high elevation. These conspire to 
produce intricate spatial variability that is hard for models to capture. A typical approach 
is to partition the space into homogeneous regions and model each region separately. While 
a number of statistical nonstationary spatial models have been proposed for univariate fields 
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Figure 1: Map of elevations in Colorado (in meters) and the 145 locations used from the 
Global Historical Climatology Network. Four locations we later use for cross-validation are 
denoted a (Kit Carson), b (Estes Park), c (Buena Vista) and d (Delta). 



(Fuentes. 


2002; 


Haas 


1990 


Hiadon. 


1998: 


Kim et al. 




2005: 


Paciorek and Schervish. 


2006; 


Pintore and Holmes 


■ 


2006 


:\sampson and Guttorp. 


1992 


Stroud et al. 


2001 


) , fewer are avail- 



able for multivariate spatial simulation which is of key concern 



and maximum temperature simulatio n ft G elf and et al. 



2012 



Shaddick and Wakefield . 120021 ) 



2004 



Jura . 



or simultaneous minimum 



2011 



Kleiber and Nychka 



Some literature in geography and the atmospheric sciences is concerned with determin- 
istic interpolation of observed weather variables, often over domains with complex ter- 



rain ( 



Daly et al 



994; 



Hiimans et al 



2005 



Hutchinson. 



1995 



Legates and Willmott 



1990 



Price et al 



2000 



Running et al. 



1987 



Thornton et al 



1997 



Willmott and Matsuurax . 



19951 ) 



The common theme among these approaches is the inclusion of high resolution digital eleva- 
tion maps as well as other physical information such as slope and aspect to deterministically 
interpolate meteorological variables. Wh ile most of these mode ls are sophisticated physical 
interpolation schemes, they (apart from \Thornton et all . 119971 ) are chiefly concerned with 



monthly or annual average quantities, and do not produce stochastic realizations of daily 
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weather, which is our primary interest. These schemes are also typically ad hoc, and are not 
based on a formal statistical model. 

Stochastic interpolation and simulation of physical variables has persistent interest in the 
statistics literature. Often, precipitation holds the prim ary interest as its mixed discrete- 



continuous and skewed nature pose substantia 



2003 



Brown et al. 



2001 



Durban and Glasbey 



challenges hAilliot et al 



2001 



2009 



Hughes and Guttorp 



Allcroft and Glasbew. 



1999 



Sanso and Guenni 



2000). However, rec ent authors have acknow 



edged the difficulties of tempera ture modeling 



2006|), and 



Gelfand et al\ ( 120051 ) is one of few 



in complex terrain \Paciorek and Schervish . 
to simultaneously model temperature and precipitation. 

The study domain in this paper is the state of Colorado. Figure [1] illustrates the challeng- 
ing terrain of Colorado, with eastern plains dipping to a minimum elevation of approximately 
1000m, and the Rocky Mountains of central Colorado peaking out at above 4000m. The front 
range, the ridge separating the Rocky Mountains from the eastern plains (running north- 
south on approximately the —105° longitude line) is especially difficult to accommodate using 
the currently available multivariate covariance models, most of which are isotropic models, 
and do not allow for sudden boundaries or even gradually evolving spatial structures across 
a domain. The 145 locations shown in Figure [1] are a subset of station s from the Global 



Peterson and Vose 



19971) • Daily ob- 



Historical Climatology Network Database (GHCND; 
servations of minimum and maximum temperature are available between a time period of at 
most 1893 through 2011. Associated with each observation is a quality flag provided by the 
GHCND; we removed all flagged observations to avoid poor quality observations. 

In this paper we propose a framework for bivariate stochastic temperature simulation that 
splits the model into two components. The first component represents local climate, allowing 
the average behavior of minimum and maximum temperature to vary with location, which 
is of critical concern in regions such as Colorado with the average behavior of temperature 
in the Rocky Mountains being vastly different than that over the eastern plains. The second 
component can be interpreted as daily weather, yielding local variability in space and time, 
and preserving the spatial correlation between both processes. 
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2 Stochastic Model 



Consider the bivariate process of minimum temperature, Z^(s,t), and maximum tempera- 
ture, Zx(s,t), at location s e R 2 on day t = 1, . . . ,T. Our model for the bivariate process 
is 

Z N (8,t)=P N (8)'X N (8,t)+W N (8,t) (1) 

Z x (s, t) = /3 x ( S yX x (s, t) + W x (s, t). (2) 

The vector of coefficients fl^s) = (/3 0i (s), (3u(s), . . . , (3 pi (s)Y, for i — N,X, may be of differ- 
ent length for minimum and maximum temperature, allowing for distinct sets of covariates, 
although for notational simplicity we assume both processes share the same number of covari- 
ates, p+1. The covariates Xj(s, t) = (X 0i (s, t), . . . , X pi (s, t))' typically involve autoregressive 
and seasonality terms, and if available, can contain additional information such as regional 
climate model output. It is convenient to view the model of ([T]) and ([2]) as a sum of "local 
climate" plus "weather" . The local climate is dependent on spatially and temporally varying 
covariates, and whose coefficients vary across the domain allowing for the relative influence 
of each covariate to depend on location. The weather terms, W^(s, t) and W x (s, t), capture 
small scale variability and correlate the bivariate temperature process across space. 



2.1 Local Climate Component 



The coefficients Pki{s), for i = N,X and k — 0, . . . , p, allow the average behavior of tempera- 
ture to vary with location. This is crucially important in areas of complex terrain or over large 
domai ns where variable orography and general circulation patt e rns give rise to varying cli- 



mate ( Chandler 



2005;\Johnson et al. 



2000 



Kleiber et al. 



2012h . \Pepin and LoslebeA fj2Q02h 



point out that climate change trends in Colorado are highly dependent on the terrain. Direct 
estimates of these coefficients are usually only available at locations within the observation 
network, so we model the coefficients as spatial Gaussian processes. In particular, we suppose 
(3ki(s) has mean fiki and Matern covariance augmented with a nugget effect, with varianc e 
parameter range a^i, smoothness and nugget effect \Guttorp and Gne^mql . 120061 ). 
The goal of a spatial model for the coefficients (3ki{s) is for interpolation from the observa- 
tional network locations to a chosen grid. The Matern is an isotropic covariance function 
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that is especially useful for kriging \Steim , Il999l ). One might consider using a nonstationary 
function for the coefficient covariance model, but in our experience (see the example section 
below), the simpler stationary model works well for local climate interpolation. 

For Colorado, we use the following covariates: 



Xjv(s, £)=(!, cos 



2ttA 
365 J 



sin 



2nt 
365 



Z x (s,t-l),Z N {s,t-l),r t 



(3) 



with the corresponding case for Xj(s, t) reversing indices N and X. The harmonics allow for 
seasonality in minimum and maximum temperature, and we include bivariate autoregressive 
terms to account for temporal persistence of temperature. The final covariate, r t is a linear 
drift of length T between —1 and 1 (for numerical stability), which we include to control for 
temperature trends over the 119 year period of our dataset, noting that these trends do not 
necessarily reflect global warming. These covariates were selected using a BIC criterion at all 
individual stations; that is, fitting a model to each location independently, the model with 
all of the above covariates had the smallest BIC value for all stations within the GHCND 
in Colorado, as compared to any subset of the selected covariates. We considered models 
with higher order harmonics and autoregressive lags, but the results were nearly identical to 
those presented below, hence we favor the simpler set of covariates. 

Suppose we observe the bivariate process (Zn(s, t), Zx{s, t))' at locations s — a\, . . . , s n 
and time points t = t\, . . . , tx- At each location within the observation network, we estimate 
local parameters (3ki(s) by ordinary least squares. These estimates have low uncertainty, 
as in the Colorado network, the location with the sparsest observational record still has 
more than 10,000 available observations. Conditional on the estimates ^(s), we estimate 
the spatial Gaussian process parameters fiki, clm, ^k% and by maximum likelihood, 
exploiting the Gauss ian process assumption. These spatially varying coefficients models 

rave been used for p robabilistic forecasting, with a similar two-step 



\Gelfand et all 1200311 
estimation procedure uKleiber et al. 



2011a.b 



At an arbitrary location s , not necessarily within the observa tion network, we spatially 
interpolate the estimates f3 ki = 0ki{ 3 i)i ■ ■ ■ > fiki{ s n))' via kriging ( Cressiei Il993 ). In partic- 
ular, the kriging estimator is 
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and the interpolation variance is 



a 



where c! = (Cov(/3 fci (s ), Pki(si)), Cov(/3 fci (s ), 0ki(s n ))) and (S)^ = Cov(/3 fei (s,), p k i(s e )) 
for j,£ = 1, . . . ,n. As kriging is an exact interpolator, when so — s ^ for any i = 1, . . . , n, 
the interpolator returns the ordinary least squares estimate (3ki{s£). 

In the next section, we exploit a nonparametric estimator of the covariance function for 
the bivariate weather process. Key to t he nonparametric estimator being consistent is a large 
number of realizations of the process uKleiber and Nychkal |2012| ) which we have available 
for the residual weather processes, whereas the coefficient processes of the local climate 
component have only one realization. Hence, we favor the parametric model with two-step 
estimation procedure for local climate. 



2.2 Weather Component 

To simulate spatially correlated fields of minimum and maximum temperature consistent 
with observed spatial patterns, we require a bivariate spatial model for Wn(s, t) and Wx{s, t). 
In particular, we model these weather processes as a zero-mean bivariate spatial Gaussian 
process indexed by day of the year. For locations x, y, and arbitrary time point t, the 
bivariate covariance model is 

Cov(W i (x,t),W j {y,t + l)) = Q (4) 
Cov(Wi(x, t), W t {y, t)) = Cu(x, y, d(t)) + n{x, y) 2 l [x=y] (5) 
Cov(W t (x,t),W,(y,t)) = C tl (x,y,d(t)) for i^j (6) 

for i, j = N, X where d(t) G {1, . . . , 365} is just the calendar day of the year on which time 
point t falls. 

The covariance model of fl4]), (jSJ) and ([6]) implies some important assumptions. First, we 
assume temporal dependence has been accounted for in the local mean function (for instance 
via autoregressive terms) so that the weather process is temporally independent, hence fl4]). 
Indeed, exploratory plots such as autocorrelation functions and empirical covariance func- 
tions indicate the bivariate autoregression of (j3J) is sufficient to account for the temporal 

7 



persistence of temperature in Colorado, see the example section below. Second, the covari- 
ance and cross- covariance functions Cu(x,y,d(t)) and Cij(x,y,d(t)) depend on the day of 
year, allowing the bivariate process to have seasonally dependent second-order structure. 

In (E]), Ti(x) 2 = Ti(x,x) 2 is a local nugget effect, accounting for small scale variability 
as well as measurement error. In the geostatistical literature, Cu(x,x,d(t)) is often termed 
the marginal varian ce, while Cu(x , x, d(t)) + Ti(x) 2 is called the sill, i.e. , the total variance 



Christensen 



2011 



at a given location ( 1 Cressid Il993l ). Unlike most geostatistical models ( 
being a notable departure), we allow the nugget effect to vary with location, as we expect 
the small scale variability to be highly dependent on orography. 

At any fixed time point t (i.e. calendar day d(t)), we require the matrix- valued covariance 
function 



C(x,y,d(t)) 



C NN {x,y,d(t)) C NX (x,y,d(t)) 
C XN (x,y,d(t)) C xx (x,y,d(t)) 



(7) 



to be a nonnegative definite matrix function. Specifically, at arbitrary locations s 1; . . . , s n , 
the covariance matrix of the random vector 

(W N ( Sl , t), W x (s u t), W N (s 2 , t), W x (s 2 , t), . . . , W N (s n , t), W x (s n , £))', 

which is made up of blocks C(s&, S£, d(t)), must be nonnegative definite. 

Over regio ns with complex terrain, temp erature observations can exhibit substantial non- 



stationarity uPaciorek and Schervish 



20061 ) . Wh ile some multivariate spatial models that 



can a ccount for nonstationarity are available (e.g. \ Gelfand et al 



2004 



Kleiber and Nychka 



20121 ). these are parametric models with locally varying parameter functions that are dif- 
ficult to estimate. We aim to exploit the large number of replications and reasonably well 
covered observation network of the GHCND over Colorado, and propose a nonparametric 
estimator of the matrix-valued covariance function that retains nonnegative definiteness. In 
particular, suppose the bivariate process is observed at locations s^, k = 1, . . . ,n and times 
t — 1, . . . , T. Then our nonparametric estimator of Cij(x, y, d(to)) in (JSJ) and (jSJ), at location 
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pair (x, y) and time point to is 

C ij (x,y,d(t )) = 

EL ELi ELi K Xt (IM(tp), d(t)\\ d ) K x (\\x-s k \\)K x (\\y- s e \\) W t (s k , t)W 3 {s^ t) 
EL ELi ELi ^(||d(to), d(*)||d)ifA(||a! - s.lD^dlj, - S ,||) 

(8) 

for i,j = N,X. Here, K x is a kernel function with bandwidth A, we use i^dl^H) — 
(1/A) exp(— ||h,||/A). We use the Euclidean norm || • ||, and the distance function is 
the distance between days of the year so that ||g?i, d 2 ||<2 = \d\ — d 2 \ for \di — d 2 \ < 182 
and 1 1 c?i , cZ 2 ] I = 1 365 — \d\ — d 2 \\ for \d\ — d 2 \ > 182, where di,d 2 = 1, ...,365, for ex- 
ample 1 1 1, 365 (I = 1. Occasionally Zi(s k ,t) (and subsequently Wi(s k ,t)) is not available 
in practice due to instrument failure or disruptions in communications. The estimator 
we use operationally is a slightly modified version of OH]), where we make the convention 
Wi(s, t)tr w .r Sit \ is observed] = when Wi(s,t) is missing. It is convenient to define the single- 
time-point smoothed empirical covariance function 

Rij(x,y,t) = 

V — s i\\) Wi(Sk, t)Wj[S£, t)t[w t (s k ,t) is observed] ~&-[Wj(s e ,t) is observed] 
Efc=l J2e=lK\(\\ X ~ s k\\)K X (\\y - S^Dl^s^is observed] l[W 3 (s«,t) is observed] 

(9) 

Notice Rij(x, y, t) is just a (spatially) smoothed empirical covariance function over the avail- 
able observations on day t. Rij(x,y,t) is a nonnegative definite multivariate covariance 
function, a property we show in the Appendix. Our adjusted version of (JE]) that accounts 
for missing observations then is 

Eti Kx t (\\d(t ),d(t)\\ d ) R^ x, y, t) 
EliK Xt (\\d(t ),d(t)\ 

As Cij(x,y,d(to)) is a positively weighted linear combination of multivariate covariance 
functions, it is again nonnegative definite. The estimator for missing observations ( flOl) 
reduces to the original estimator ([8]) when no observations are missing, and hence ([8]) is also 
nonnegative definite. 

The estimator (jSJ) is a smoothed version of daily empirical covariance matrices. The 
first level of smoothing yields an estimate of spatial covariance at any arbitrary location 



c^x, y, d(to)) = ^ j^^rr/jTx ,;.r„ r • ( iq ) 
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pairs in the domain. The temporal smoothing shares information between adjacent time 
points, where we assume that spatial covariance on a given day is similar to that in a short 
period leading up to that day, and in a short period following that day. This esti mator is 
a gene r alization of k e rnel s moo thed empirical c ovariance estimators considered by 



(1993). 



Guillot et al. 



fl200lh and 



Oehlert 



Jun et al 



(120 111 ) to the multivariate process setting evolving 



across time. 



We estimate the time bandwidth \ t by predictive leave-one-out cross-validation, leaving 
out local empirical variance estimates. The estimated bandwidth for time is \ t = 7.8 days. 
We use cross-validation for the temporal bandwidth as we assume the temporal evolution of 
spatial covariance is slowly evolving across time, e.g. we do not expect a sharp change in spa- 
tial covariance between June 1 and June 2. In our experience, using cross-validation for the 
spatial bandwidth parameter A oversmooths the spatial covariance function. When kernel 
smoothing a mean function, cross-validation is generally acknowledged to yield more vari- 
ability than is ex pected for a smoo t hly v arying mean function, and typically the bandwidth 



must be inflated (I Wand and Jones 



1995|) . However, this experience is under the assumption 



that the mean function is varying smoothly across the domain, and in regions of complex 
terrain we expect the opposite behavior, where sharp boundaries of the covariance function 
may exist due to sudden changes in elevation. For example, cross-validation implies the 
optimal spatial bandwidth is 75 km, which implies an effective range of the kernel function 
(that is, up to 5% weight) of approximately 225 km, greatly oversmoothing regions such as 
the San Luis Valley in southern Colorado, at approximately 100 km across. Hence, we choose 
a bandwidth such that the effective distance of the kernel function coincides with the 5% 
quantile of all intersite distances (62 km); the heuristic argument is that, for approximately 
evenly distributed observation locations, the covariance estimator at a given location uses 
the nearest 5% of available network locations, and down-weights remote locations; this ad 
hoc criterion implies a spatial bandwidth of A = 22 km. 

The estimator ([8]) is asymptotically unbiased for Cij(x, y, d(to)) when the domain sample 
size increases and the bandwidth decreases to zero sufficiently quickly. A short argument 
is given in the Appendix. In fact, it can be sho wn that the estimator is co nsistent for 
Cij(x,y,d(t )), using arguments similar to those o uKleiber and Nychkal (120121 ). but this is 
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beyond the scope of the present paper. 

All that remains to be estimated is the local nugget effect Tj(s) 2 . At each observation 
location Sk,k = 1, . . . , n, and time point t — 1, . . . , T, let Wi(sk, t) be the estimated residual 
t) — / 9 i (s ( ! C ) / Xj(s/ c , t). Define the local empirical variance on day d — 1, . . . , 365 as 

diM2 = mm-Q E w ' {Sk,tf 

7ri 1 w J {t\d(t)=d} 

where # denotes cardinality of the set, with the natural redefinition for missing values of 
Wi(s k ,t). Intuitively, a good estimator for Ti(s k ) 2 is 

^ 365 

Ti(s k ) 2 = — (d-i(s k , d) 2 - Cu(s k , s k , d)^ (11) 

d=l 

since, by the law of large numbers, &i(s k , d) 2 — > Ca(s k , s k , d) +Ti(s k ) 2 where the convergence 
is taken as T — > oo, and by the argument in the Appendix, Cu(s k , s k ,d) — > Cu(s k , s k ,d). 
While theoretically appealing, in practice due to the smoothing in Cu, at some locations the 
estimate fi(s k ) is negative. Hence, in similar spirit we use f[TT|) . but set the invalid estimates 
to zero. 

Estimates of Tj(s) 2 are gathered at arbitrary locations, i.e., not necessarily within the 
observation network, by imposing a probabilistic spatial structure on Tj(s). In particular, 
we model Tj(s) as a Gaussian process with spatially constant mean and Matern covariance 
function, augmented with a nugget effect. Just as for the spatial parameters of the (3 ki (s), we 
estimate the spatial parameters of Tj(s) by maximum likelihood, conditional on the estimates 
{^i( s k)}k=v While the estimates fi(s k ) at observation locations are always valid, the kriging 
interpolator of Tj(s) may occasionally take on very small negative values; in our example 
below we did not experience such an issue, but in other domains these degenerate estimates 
may be artificially set to zero. 



3 Minimum and Maximum Temperature in Colorado 

We fit our model to the data from the 145 GHCND locations shown in Figure [TJ For 
simplicity, we removed all leap days from the 119 years of available data, so that each year has 
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Figure 2: Empirical autocorrelation functions for minimum and maximum temperature resid- 
uals at locations (a) Kit Carson, (b) Estes Park, (c) Buena Vista and (d) Delta. 



365 days. Using all available data, we fit local climate parameters by ordinary least squares, 
and estimate temporally varying multivariate spatial covariances using the nonparametric 
estimator (TTUT) applied to the observed residuals. We then simulate the bivariate process for 
a 119 year trajectory to compare to the observed bivariate series. The first day's (January 
1, 1893) simulation requires autoregressive terms in (OQ) and ([2]); we initialize using the 
climatological domain average of minimum and maximum temperature on December 31. The 
resulting simulations are masked to share the same missing value pattern as the observations. 

Recall the assumption implied by Equation (jl]), where we assume temporal dependence 
has been accounted for in the local mean function via the bivariate autoregression. Figured 
contains empirical autocorrelation functions for the observed residuals Wn{s, t) and Wx{s, t) 
at four network stations, shown in Figure [TJ These locations we view as representative of 
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Figure 3: Scatterplots comparing empirical pairwise station correlation to simulated corre- 
lations using a bivariate stationary model (grey dots) or the nonstationary nonparametric 
model (black dots) over the summer months (JJA). The diagonal line indicates perfect agree- 
ment between model and empirical correlations. 

four distinct regimes of Colorado: eastern plains (a, Kit Carson), front range (b, Estes Park), 
Rocky Mountains (c, Buena Vista) and the western slopes (d, Delta). It is evident that the 
bivariate autoregression accounts for the majority of temporal persistence in temperature; 
the maximal lag-1 autocorrelation coefficient for the residual processes at these four stations 
is 0.06 at Estes Park, whereas all other coefficients are less than or equal to 0.03. 

To motivate the flexibility of the nonparametric estimator (JSJ), we compare it to a state- 



of-the-art isotropic bivariate spatial mod e 



nGneitina et al 



2010 



Apanasovich et al 



. In particular, we fit a bivariate Matern model 



20121 ) augmented with a nugget effect, where 



Ca(x,y,t) 
C NX (x,y,t) 



2^- 1 r(z/,- 



a,- as 



y\\TK Vi {ai\\x 



Pnx: 



y\\) + rtl [x=y] 
a N x\\x - y\\Y l K l/NX (a NX \\x - y\ 



(12) 
(13) 



2-nx~it(u nx ) 

for % = N,X, where K u is a modified Bessel function of the second kind of order u, and 
v nx = {y N + z/x)/2 and a X x = min(ajv, ax)- We fit the parameters by maximum likelihood, 
viewing each bivariate estimated residual (Wjv(s, t), Wx(s, t)) as independent across time. 
In the stochastic weather simulation literature, it is customary to fit separate models for 
each season. While our nonparametric estimator is available on any day, to facilitate com- 
parisons to the bivariate Matern, we fit both models to only the summer months (JJA), and 
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Figure 4: Plots of spatial correlation and cross-correlation functions, CV,(s , -, d(t)), on d(t) = 
June 1 where s is a grid location in the eastern plains (top row) or a grid location in the 
Rocky Mountains (bottom row), with grid locations indicated by black dots. Each pixel's 
color indicates the model estimated spatial correlation between the pixel location and the 
dot. 

compare empirical to simulated correlations and cross-correlations under both the isotropic 
and nonparametric models; Figure [3] displays these results. The stationary model tends 
to overestimate spatial correlation for both minimum and maximum temperature, whereas 
our nonstationary model adequately captures low and high correlations simultaneously. The 
third panel of Figure E] shows empirical against simulated cross-correlations. Substantial 
nonstationarity of cross-correlation across Colorado is well modeled by our nonparametric 
approach, but the stationary model clearly fails, putting most cross-correlations at around 
0.10, whereas the empirical estimates suggest the true cross-correlations should vary between 
-0.10 and 0.40. 

Our nonparametric matrix covariance estimator (jSJ) accommodates nonstationary behav- 
ior of the multivariate process. Figure H] shows two covariance functions on June 1, one whose 
first argument is based at a grid location in the eastern plains of Colorado, and the second 
covariance function whose first argument is based at a grid location in the Rocky Moun- 
tains. The top row is the covariance function for the plains based grid location; particularly 
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Figure 5: Same as Figure HI except for d(t) = January 1 instead of June 1. 
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for maximum temperature, and lesser so for minimum temperature, there is strong positive 
within variable correlation throughout the plains region, suggesting that maximum temper- 
atures are highly correlated across the plains. At the front range boundary (approximately 
—105° longitude), there is a sharp drop off in spatial correlation from approximately 0.80 
over the plains to 0.40 in the Rocky Mountains. This is due to the fact that temperature 
is more highly correlated within the two main types of topography of Colorado, either the 
plains or mountains, but not between the two types. Hence, our estimator is able to capture 
the sharp boundary between the eastern plains and Rocky Mountains for within variable 
spatial correlation. Our estimator also identifies the positive cross-correlation between mini- 
mum and maximum temperature in the plains, but allows the two processes to be effectively 
independent over the Rocky Mountains. This nonstationarity of cross-correlation is very 
difficult to acc ommodate using extant m odels, and has only been recently acknowledged in 



the literature ttKleiber and Genton 



20121 ). 



Not only does our estimator allow for substantial nonstationarity, the amount and type of 
nonstationarity is allowed to vary across time. Figure [5] shows the same plots of spatial direct 
and cross-correlation on January 1, during winter, as opposed to the summer estimates of 
Figure HI In terms of direct covariance, we see the length scale of minimum temperature 
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Figure 6: Estimated nuggets f N (s) and f x (s) at observation network locations, units are 
degrees Celsius. 

correlation drastically increase for both the plains and mountain based grid locations. In 
the plains, the spatial correlation structure of maximum temperature is similar during both 
the winter and summer; on the other hand, this spatial correlation in the mountainous 
region over winter has a substantially different pattern than over summer. The correlation 
structure of the weather component for maximum temperature in the Rocky Mountains 
is clearly nonstationary, implying lower correlation between the example grid point and 
the southwestern slopes of the Rockies, but having higher correlation along a northwest to 
southeast transect along the western slopes and through the Rocky Mountains; this pattern 
makes sense climatologically, as the band of high correlation connects the low lying western 
Grand Valley area through the lower mountains north of the San Juan chain to the San Luis 
Valley in southern Colorado. A similar pattern is present for the cross-correlation function, 
which is distinct from the summer behavior which indicated near-independence between 
minimum and maximum temperature over the complex topography. 

A notable departure of our model from typical geostatistical approaches is in allowing the 
nugget effect to vary with location. Our motivation is that the small scale spatial structure is 
expected to be dampened in the eastern plains with stable orography, but potentially inflated 
over the mountainous region of Colorado. Figure [6] displays the local estimates fj(s) for 
i = N, X at the locations within our observation network. For both minimum and maximum 
temperature, the nugget effects tend to be less over the eastern plains, indicating less fine 
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scale spatial structure (although there is yet some evidence of small scale structure in the 
maximum temperature nuggets here). Over the Rocky Mountains, especially the northern 
Rockies, minimum temperature exhibits inflated nugget effects, indicating fine scale spatial 
processes in the complex terrain. Similarly, the finest scale spatial structures indicated by 
these nugget effects for maximum temperature falls almost directly along the front range, the 
longitude line of approximately —105°, indicating highly variable maximum temperatures 
between the boundary of the plains and sudden mountainous terrain. The inclusion of a 
spatially varying Tj(s) allows the statistical model to retain increased variability along the 
front range, for example, while simultaneously generating tempered fields over the eastern 
plains, and fields of medium variability over the main Rockies and western slopes. 

An increas ingly important c onsid eration in climate science is the effect of climate change 



on extremes (\Easterlina et all 120001 ). Our model is not explicitly designed to replicate ex- 



treme events, as we focus mainly on the first and second order properties of minimum and 
maximum temperature. Figure [7] shows Q-Q plots for daily domain wide extrema. In partic- 
ular, we find the minimal and maximal domain-wide temperatures Z itUlin (t) = min s {Zj(s, £)} 
and Zj jmax (i) = max s {Zj(s, £)}, and compare simulated to observed daily statistics for 
i = N,X. Our model replicates the statistical properties of Z N!rain (t), Zx, m m{t) and Z Njraax (t) 
very well, at even the most extreme tails of these domain extrema. However, we simulate 
domain wide maximal maximum temperatures that are slightly too high, on average about 
2°C. Overall, even though our approach does not explicitly model extreme temperatures, we 
are able to capture the spatial extrema with reasonable accuracy. 

While our model adequately replicates domain wide extrema, the related quantity of 
spatially consistent local extrema is critically important to replicate. In particular, for energy 
use forecasting and modeling, if a large number of locations experience unusually low or high 
temperatures simultaneously, then the load on the energy grid can be much greater than if 
the temperature anomaly were highly localized. Figure [8] shows log frequencies (i.e. the 
log number of days) of the number of stations whose local weather process Wi(s,t) either 
exceeded the local 90% quantile (that is, the quantile using only data from location s) or fell 
below the local 10% quantile, corresponding to local hot or cold events, respectively. Our 
approach captures the spatial frequencies of unusual local cold temperatures extremely well, 
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Figure 7: Q-Q plots for daily spatial extrema, comparing (a) domain minimum of minimum 
temperature, (b) domain maximum of maximum temperature, (c) domain maximum of min- 
imum temperature and (d) domain minimum of maximum temperature, units are degrees 
Celsius. 

and tends to simulate local heat events over slightly inflated regions when many stations 
experience hot events, although usually fewer than seven extra days on average. 

Our nonparametric weather component covariance estimator (|S]) is not optimized for 
cross-validation. To assess the interpolative properties of our estimator, we hold out data 
from the four network stations shown in Figure (TJ representing four distinct regimes of 
Colorado. We predict the local standard deviations Cu(s, s, d) 1 ^ 2 for i = N,X, and compare 
these to the locally estimated values of Ca(s, s, d) 1 ^ 2 when station data is retained. Figure 
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Figure 8: Log frequency of observed and simulated local residual threshold exceedances. 
Each bar's height is the log freqeuency (i.e. log number of days) that an exact number of 
the observation network stations had weather that exceeded the local 90% quantile for (a) 
maximum temperature, or whose weather fell below the local 10% quantile for (b) minimum 
temperature. 

contains the local and predicted estimates for all days of the calendar year. Clearly the 
weather component variability is highly dependent on season as well as location, particularly 
for maximum temperature there is substantially greater variability in the eastern plains (3°- 
5°C), compared to the mountain regions (2°-3.5°C). Our predictive local standard deviations 
(dashed lines in Figure H]) generally agree closely with the local estimates, although there is 
a slightly tendency to under-predict local standard deviation at Kit Carson by 0.1°-0.3°C. 
Not only are the raw values well predicted, but the climatological curvature is preserved as 
well; for example, we successfully replicate the increased variability of maximum temperature 
over the western slopes during springtime with relatively constant variability throughout the 
three remaining seasons (panel d) while simultaneously producing significant seasonality over 
the eastern plains, with low variability during summer and high variability during winter 
(panel a). 

Table CD shows the interpolated coefficients with predictive standard deviation, along with 
the locally estimated parameters Pki{ s ) fc> r k = 0,...,5 and i = N,X for s being one 
of the four held out network stations. All locally estimated parameters are within the 
95% predictive confidence interval, except for four cases for maximum temperature. Our 
predictive intervals are calibrated; the coverage of the 95% interpolation intervals for leave- 
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Figure 9: Locally estimated standard deviations (Cu(s, s, d) 1 ^ 2 ) for i = N, X on all days of 
the calendar year d — 1, . . . ,365, and predicted standard deviations for the four hold out 
stations s = (a) Kit Carson, (b) Estes Park, (c) Buena Vista and (d) Delta. 

one-station-out cross-validation over all locations was, at worst, 92.4% for /?2x(s)- Notice 
that the local estimates vary substantially between locations, indicating that indeed the 
local climate varies over the domain. Hence we are able to successfully predict the local 
weather component parameters and local climate component parameters at these four hold 
out locations which are representative of four regimes in Colorado. 

Finally, we illustrate the final product of our approach in Figure [10] which displays 
four days of gridded simulations of minimum and maximum temperature over Colorado. 
Marginally, we visually see the temporal persistence of temperature over a period of days, 
as both minimum and maximum temperature experience a period of cooling over June 1-4. 
Notice the effect of local climate is to keep the Rocky Mountain region cooler for both vari- 
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Table 1: Interpolated estimates (with predictive standard deviation) of the local climate 
component coefficients with the validating locally estimated parameters. Locations are s = 
(a) Kit Carson, (b) Estes Park, (c) Buena Vista and (d) Delta. Predictions are starred if 
the truth is outside of the predictive 95% confidence interval. Units are degrees Celsius for 
Po, (3\ and /3 2 , unitless for (3 3 and /3 4 , and degrees Celsius per century for (3 5 . 
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ables, allowing higher minimum and maximum temperatures to fall over the eastern plains 
of Colorado. We also see slightly warmer temperatures on the western slopes, as the Rocky 
Mountains decay in elevation to the western border of Colorado. The cross-correlation be- 
tween the two variables is also present, as both variables are seen to cool across the domain 
simultaneously. 



4 Discussion 

In this paper, we introduce a framework for stochastic bivariate minimum and maximum 
temperature simulation over complex domains. The framework distinguishes between local 
climate and weather processes. The local climate is accommodated through a linear model 
whose coefficients are spatially varying, and the weather process is modeled as a bivariate 
spatial Gaussian process with a nonparametric estimate of the matrix valued covariance 
function that retains nonnegative definiteness at arbitrary locations. We successfully capture 
the temporally varying spatial dependence between minimum and maximum temperatures 
over the state of Colorado, which exhibits challenging complex terrain that is difficult for 
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Figure 10: Gridded simulation of daily minimum and maximum temperature on days June 
1-4. 



extant models to accommodate. 

Our nonparametric estimator smooths multivariate spatial covariance over space as well 
as time. This approach allows spatial dependence to be highly different during winter than 
during summer, for instance, and also retains nonstationary spatial structures both within 
each process and between processes. The estimator is available at any location, not only 
those within the observation network, and always retains nonnegative definiteness, allowing 
for gridded simulations. The estimator relies on kernel smoothed empirical covariance func- 
tions, and our current approach to spatial bandwidth selection is ad hoc. One future route 
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of research may be to decide on a quantitative approach to bandwidth selection when sharp 
boundaries and highly variable covariances are expecte d across the study doma in, notably 
different than most mean function smoothing literature ft Wand and Jonesl Il995l ). A second 
potential direction of research may be to develop a nonparametric kernel smoothed esti- 
mate of the multivariate covariance function that is robust against outliers and still retains 
nonnegative definiteness. 

While our approach does not explicitly model extremes, our simulations indicate reason- 
able replication of tail behavior, even domain-wide extrema. A limitation in using Gaussian 
processe s is th at there is a lack of clustering at high levels, both spatially and temporally 



nSibuyal 119601 ) . This is one potential explanation for the behavior of Figure [7] (b), where 



domain-wide maximal maximum temperatures were simulated slightly above the observed 
extremes, although we would expect to see similar behavior in panels (a), (c) and (d). An 
approach that includes a Gaussian process model for the bulk of the distribution along with 
a model for spatial extremes may improve extremal performance. 

One consideration of our model is that we do not explicitly force the simulation of max- 
imum temperature to be greater than or equal to minimum temperature; in our Colorado 
example maximum temperature was less than minimum temperature for approximately 
one tenth of a percent of our simulations. 



t may be of interest to adopt the models of 
Jolliffe and Hopa (119961 ) or \Jones and Larsen (12004 ) to our situation if this issue is of criti- 



cal concern. 

The clearest route of future research is to extend our ideas to a full stochastic weather sim- 
ulator that can simulate spatially correlated fields of multiple variables such as minimum and 
maximum temperature, precipitation amount, solar radiation, wind direction/wind speed 
and relative humidity simultaneously. Indeed, in complex terrains the practitioner will need 
to rely on highly flexible spatial models to replicate the strong nonstationarities exhibited 
by these various processes, as well as the complicated spatially evolving relationship between 
them. 



23 



Appendix 



In this Appendix we show the nonparametric estimator is nonnegative definite, from 
which it follows that (JSJ) and ( TTUj) are also nonnegative definite. Below, we present an 
argument that (JHD is asymptotically unbiased for Cij(x,y,d(t )). 

The nonnegative definiteness property is not restricted to a bivariate process, so as- 
sume there are p spatial processes Wi(s,t),i = 1, . . . ,p, with observation network locations 
s m ,m = l,...,n. Define Ui(s,t) = Wi(s, *)l[wi(«,t) is observed], noting that if Wi(s,t) is un- 
available at a particular location and time, Ui(s,t) = 0. 

Consider evaluating Rij(xk,xt,t) at any arbitrary locations Xj, and x#,k,£ = 1,...,N, 
and define the arbitrary vector a = (a n , . . . , a 1A r, a 2 i, . . . , a p ^). Set £ to be the covariance 
matrix made up of the functions Rij{-, -,t) corresponding to the random vector 

(Wi(a;i, *), W 1 (x 2 , t),..., W p (x N , t))'. 

Then, absorbing the denominator into the kernel functions of Rij(xk,X£,t), and writing 
R' i: j(xk,Xi,t) for this normalized function, we have 

V N 

i,j=l k,£=l 
p N n 

= ^ ^CLikCtji K *(\\ x t~ a m\\)K\(\\&k - *r\\)Ui(8 r ,t)Uj(8 m ,t) 

i,j=l k,£=l m,r=l 
n p N 

= S S K^dl 35 * - s r \\)Ui(s r ,t)) (a je K x (\\x e - s m \\)Uj(s m ,t)) 

m,r=l i,j=l k,l=l 

(n p N \ 2 

53 535ZaifciirA(||ajjfc-8 r ||)t/i(a r ,t)] > 0. 
r=l i=l fc=l / 

To show that (jSJ is asymptotically unbiased for C{j(x, y, d(t )), we disregard the smooth- 
ing over time, since asymptotically we do not have a finer resolution of time points (but 
for consistency we would assume an increasing number of realizations per each day of 
the year). In particular, suppose we observe the bivariate process (Wjv(sfc), Wx(sk)) for 
Si . . . , s n G T> C M. d , which are samples from a distribution with strictly positive probability 
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density / : V — > M + , with empirical cdf F n (x) = - Ylk=i M s k<x}> where the indicator func- 
tion is 1 if the inequality holds for all indices of x. The density /, with corresponding cdf 
F, allows the network density to vary across the domain. We additionally suppose n — > oo 
and A — > such that A ~ n~ l ' d+e for some small < e < 1/d 2 . 

Suppressing the time indexing from our notation, we can write 

n n 
k=l 1=1 

where the denominator of ([8]) is absorbed into the kernel functions of the numerator, yielding 
standardized functions K' x . Here we only consider the direct covariance estimators at a 
location x e T>\&D, Cu(x, x); the same argument applies for the direct and cross-covariance 
functions Cij(x, y) for x ^ y. We have 



EC u (x,x) 



n 2 X 



1 n 



(||aj - s k \\)K' x (\\x - 8i \\) M(Wi(8 k )Wi(8 t )) 



k,e=i 

n 



n 2 X 



+ 



\&Y l K' x {\\x-8 k \\)K' x {\\x- aSC^au, s t ) 



kl=l 



n 2 X 



n 



\K\\x-a k \ 



n(s k , s k ) 2 . 



(14) 
(15) 
(16) 



_t=L 



Invoking Lemma 7 oi xKleiber and Nychkm (120121 ) . in the limit as n — > oo, we can pass from 
the sum to the integral. Assume the empirical cdf F n is close to the limiting cdf F, where 
sup x \F n (x)—F(x)\ = D n where D n = o{l/{nX d )). This rate holds, for example, if V — [0, 1], 
F is the uniform density and F n is the empirical cdf of the uniform grid (l/n,2/n, . . . ,n/n). 
For d > 1, if n grows as m d , a rate of D n ~ l/n 1 ^ can be derived for sampl ing locations on 



a regular grid with limiting uniform distribution \Kleiber and A^i/c^fcal . l2012h . Then we have 

1 



ECu(x,x) 



X 2d 
+ 



K' x (|| U - x\\) K' x (||v - x\\) C u (u, v)dFHdF(v) 



T> 2 



±Z [ K' x {\\u-x\\fr^ufdF{u) + 0{D n ). 

Making the change of variables to a = (u — x)/X and b = (v — x)/X yields 

l\ K' (||o||) K' (||b||) CuiXa + x, Xb + a;)dF(a)dF(b) 
J Jv 2 



(17) 



+ 



J K'(\\a\\) 2 r t (Xa + x,Xa + x) 2 dF(a) + 0(D n ) 
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for an appropriate translated domain V . As A ~ n ~ ^ d+£ , the second term of ([TBI converges 
to zero. The arguments from \Kleiber and Nychkax (120121 ) applied to the first term of (ITS1) 
then yield the unbiasedness of Cu(x, x). 
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